function results = olsgmm(lhv,rhv,lags,weight,inversion)

%       function olsgmm does ols regressions with gmm corrected standard errors
%       Inputs:
%       lhv T x N vector, left hand variable data
%       rhv T x K matrix, right hand variable data
%       lags number of lags to include in GMM corrected standard errors
%       weight: 1 for newey-west weighting
%               0 for even weighting
%               -1 skips standard error computations. This speeds the program up a lot; used inside monte carlos where only estimates are needed
%       inversion: determines whether the Backslash (matrix left division) operator is used or pinv(A)*B.
%                  If A is an n-by-n matrix and B is a column vector with n components, or a matrix with several such columns,
%                  then X = A\B is the solution to the equation AX = B. A warning message will be displayed if A is badly scaled
%                  or nearly singular. If A is an m-by-n matrix with m ~= n and B is a column vector with m
%                  components, or a matrix with several such columns, then X = A\B is the solution in the least squares
%                  sense to the under- or overdetermined system of equations AX = B. The effective rank, k, of A is
%                  determined from the QR decomposition with pivoting. A solution X is computed that has at most k
%                  nonzero components per column. If k < n, this is usually not the same solution as pinv(A)*B,
%                  which is the least squares solution with the smallest norm
%                  Choose inversion = 1 (default) to use the backslash
%                  operator or inversion = 2 for a pseudo inverse
%
%  NOTE: you must make one column of rhv a vector of ones if you want a constant.
%        should the covariance matrix estimate take out sample means?
%  v:           variance covariance matrix of estimated parameters. If there are many y variables, the vcv are stacked vertically
%
%  Note: program checks whether first is a constant and ignores that one for test
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------
% RETURNS: a structure
%        results.meth  = 'OLSGMM'.
%        results.beta  = regression coefficients K x 1 vector of coefficients
%        results.se    = K x 1 matrix standard errors of parameters. (Note this will be negative if variance comes out negative)
%        results.tstat = GMM tstats
%        results.F     = F TEST [Chi squared statistic    degrees of freedom     pvalue] for all coeffs jointly zero.
%        results.yhat  = yhat.
%        results.resid = residuals.
%        results.rsqr  = rsquared.
%        results.rbar  = rbar-squared.
%        results.nobs  = nobs.
%        results.nvar  = nvars.
%        results.s2  = Error Variance
%        results.dw    = Durbin-Watson Statistic
% --------------------------------------------------


if nargin == 4
    inversion = 1;
end

if size(rhv,1) ~= size(lhv,1);
    disp('olsgmm: left and right sides must have same number of rows. Current rows are');
    size(lhv)
    size(rhv)
end;

T = size(lhv,1);
N = size(lhv,2);
K = size(rhv,2);
sebv = zeros(K,N);
Exxprim = pinv((rhv'*rhv)/T);

% decide how to project left on right hand side:
if inversion == 1
    bv = rhv\lhv;
elseif inversion == 2
    bv = pinv(rhv)*lhv;
end

if weight == -1;  % skip se's if you don't want them.  returns something so won't get error message
    sebv=NaN;
    R2v=NaN;
    R2vadj=NaN;
    v=NaN;
    F=NaN;
else
    errv = lhv-rhv*bv;                  % compute residual / if there are multiple LHVs then this residuals are in columns
    s2 = mean(errv.^2);                 % compute variance of residual /
    vary = lhv - ones(T,1)*mean(lhv);   % compute sample variance of observations
    vary = mean(vary.^2);               % ''
    
    R2v = (1-s2./vary)';                                    % compute R^2
    R2vadj= (1 - (s2./vary)*(T-1)/(T-K))';       % compute adjusted R^2
    
    %compute GMM standard errors
    for indx = 1:N;
        err=errv(:,indx);
        inner = (rhv.*(err*ones(1,K)))'*(rhv.*(err*ones(1,K)))/T;   % compute S matrix in GMM formulae
        
        for jindx = (1:lags);
            inneradd = (rhv(1:T-jindx,:).*(err(1:T-jindx)*ones(1,K)))'...
                *(rhv(1+jindx:T,:).*(err(1+jindx:T)*ones(1,K)))/T;        % compute lags of S ...
            inner = inner + (1-weight*jindx/(lags+1))*(inneradd + inneradd');   % compute SUM_{j=-lag}^{j = +lag}
        end;
        
        
        varb = 1/T*Exxprim*inner*Exxprim;               % compute variance of beta with Newey-West correction to S
        
        % F test for all coeffs (except constant) zero -- actually chi2 test
        if rhv(:,1) == ones(size(rhv,1),1)
            
            varBETAS = varb;
            chi2val = bv(2:end,indx)'*pinv(varb(2:end,2:end))*bv(2:end,indx);            %compute Chi^2 test statistic alpha'*cov[alpha]*alpha
            dof = size(bv(2:end,1),1);                                                  % degress of freedom
            pval = 1-cdf('chi2',chi2val, dof);                                          % compute p-value
            F(indx,1:3) = [chi2val dof pval];
            
        else
            varBETAS = varb;
            chi2val = bv(:,indx)'*pinv(varb)*bv(:,indx);
            dof = size(bv(:,1),1);
            pval = 1-cdf('chi2',chi2val, dof);
            F(indx,1:3) = [chi2val dof pval];
        end;
        
        if indx == 1;
            v = varb;
        else
            v = [v; varb ];
        end;
        
        seb = diag(varb);
        seb = (abs(seb).^0.5);
        sebv(:,indx) = seb;
    end;
end;


% Compute Log-Likelihood value 
sigma2_ML = errv'*errv/T; % NB ML estimator of variance is biased but consistent. 
LL = -T/2 * (log(2*pi) + log(sigma2_ML)) - 0.50 * (1/sigma2_ML) * (errv'*errv);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results.meth    = 'olsgmm';
results.nobs    = T;
results.nvar    = K;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results.beta    = bv;
results.se      = sebv;
results.tstat    = bv./sebv;
results.F       = F;
results.vcov   = varBETAS;
results.rsqr    = R2v;
results.rbar    = R2vadj;

results.pval    = pval;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results.yhat    = rhv*bv;
results.resid   = errv;
results.s2 = errv'*errv/(T-K);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ediff = results.resid(2:T) - results.resid(1:T-1);
results.dw = (ediff'*ediff)/(errv'*errv); % durbin-watson
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results.LL = LL; % log-likelihood

end


